Problem Statement
In this capstone project, the goal is to build a pneumonia detection system, to locate the position of inflammation in an image. Specifically, your algorithm needs to automatically locate lung opacities on chest radiographs. Tissues with sparse material, such as lungs which are full of air, do not absorb the X-rays and appear black in the image. Dense tissues such as bones absorb X-rays and appear white in the image. While we are theoretically detecting “lung opacities”, there are lung opacities that are not pneumonia related. In the data, some of these are labeled “Not Normal No Lung Opacity”. This extra third class indicates that while pneumonia was determined not to be present, there was nonetheless some type of abnormality on the image and oftentimes this finding may mimic the appearance of true pneumonia.
Dicom original images
Medical images are stored in a special format called DICOM files (*.dcm). They contain a combination of header metadata as well as underlying raw image arrays for pixel data.
Here’s the backstory and why solving this problem matters
Data Source: https://www.kaggle.com/c/rsna-pneumonia-detection-challenge/data
Reference: https://www.who.int/news-room/fact-sheets/detail/pneumonia
!pip install pydicom
import numpy as np
import pydicom
import pandas as pd
import seaborn as sns
import matplotlib as mpl
import multiprocessing
from tqdm import tqdm_notebook
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
pd.set_option('display.float_format', lambda x: '%.3f' % x)
sns.set_style('whitegrid')
sns.set_context('notebook')
plt.rcParams['figure.figsize'] = [14, 8]
plt.rcParams['lines.linewidth'] = 2.5
from google.colab import drive
drive.mount('/content/drive/')
import os
path = '/content/drive/MyDrive/Colab Notebooks/GL-AIML/Capstone Project'
os.chdir(path)
!pwd
ls
# from zipfile import ZipFile
# with ZipFile('rsna-pneumonia-detection-challenge.zip','r') as z:
# z.extractall()
class_info_df = pd.read_csv('stage_2_detailed_class_info.csv')
train_labels_df = pd.read_csv('stage_2_train_labels.csv')
train_labels_df.info()
class_info_df.info()
print(f"Detailed class info - rows: {class_info_df.shape[0]}, columns: {class_info_df.shape[1]}")
print(f"Train labels - rows: {train_labels_df.shape[0]}, columns: {train_labels_df.shape[1]}")
Let's explore the two loaded files. We will take out a 5 rows samples from each dataset.
class_info_df.sample(5)
train_labels_df.sample(5)
print('Number of rows (unique boxes per patient) in main train dataset:', train_labels_df.shape[0])
print('Number of unique patient IDs:', train_labels_df['patientId'].nunique())
In class detailed info dataset are given the detailed information about the type of positive or negative class associated with a certain patient.
In train labels dataset are given the patient ID and the window (x min, y min, width and height of the bounding box containing evidence of pneumonia.)
def missing_data(data):
total = data.isnull().sum().sort_values(ascending = False)
percent = (data.isnull().sum()/data.isnull().count()*100).sort_values(ascending = False)
return np.transpose(pd.concat([total, percent], axis=1, keys=['Total', 'Percent']))
missing_data(train_labels_df)
missing_data(class_info_df)
The percent missing for x,y, height and width in train labels represents the percent of the target 0 (not Lung opacity).
Let's check the class distribution from class detailed info.
plt.rcParams['figure.figsize'] = [8, 6]
plt.rcParams['lines.linewidth'] = 1.5
def countPlot(target="Target"):
ax = sns.countplot(x=target, data=train_labels_df, palette="pastel")
for p in ax.patches:
ax.annotate('{}'.format(p.get_height()), (p.get_x()+0.375, p.get_height()+0.15), ha='center')
plt.show()
countPlot("Target")
f, ax = plt.subplots(1,1, figsize=(10,8))
total = float(len(class_info_df))
sns.countplot(class_info_df['class'],order = class_info_df['class'].value_counts().index, palette='Set3')
for p in ax.patches:
height = p.get_height()
ax.text(p.get_x()+p.get_width()/2.,
height + 6,
'{:1.2f}%'.format(100*height/total),
ha="center")
plt.show()
Let's look into more details to the classes.
def get_feature_distribution(data, feature):
# Get the count for each label
label_counts = data[feature].value_counts()
# Get total number of samples
total_samples = len(data)
# Count the number of items in each class
print("Feature: {}".format(feature))
for i in range(len(label_counts)):
label = label_counts.index[i]
count = label_counts.values[i]
percent = int((count / total_samples) * 10000) / 100
print("{:<30s}: {} or {}%".format(label, count, percent))
get_feature_distribution(class_info_df, 'class')
No Lung Opacity / Not Normal and Normal have together the same percent (69.077%) as the percent of missing values for target window in class details information.
In the train set, the percent of data with value for Target = 1 is therefore 30.92%.
Let's merge now the two datasets, using Patient ID as the merge criteria.
assert train_labels_df['patientId'].values.tolist() == class_info_df['patientId'].values.tolist(), 'PatientId columns are different.'
train_class_df = pd.concat([train_labels_df, class_info_df.drop(labels=['patientId'], axis=1)], axis=1)
train_class_df.shape
train_class_df.info()
Let's plot the number of examinations for each class detected, grouped by Target value.
fig, ax = plt.subplots(nrows=1,figsize=(12,6))
tmp = train_class_df.groupby('Target')['class'].value_counts()
df = pd.DataFrame(data={'Exams': tmp.values}, index=tmp.index).reset_index()
sns.barplot(ax=ax,x = 'Target', y='Exams',hue='class',data=df, palette='Set3')
plt.title("Chest exams class and Target")
plt.show()
All chest examinations with Target = 1 (pathology detected) associated with class: Lung Opacity.
The chest examinations with Target = 0 (no pathology detected) are either of class: Normal or class: No Lung Opacity / Not Normal.
For the class Lung Opacity, corresponding to values of Target = 1, we plot the density of x, y, width and height.
target1 = train_class_df[train_class_df['Target']==1]
sns.set_style('whitegrid')
plt.figure()
fig, ax = plt.subplots(2,2,figsize=(12,12))
sns.distplot(target1['x'],kde=True,bins=50, color="red", ax=ax[0,0])
sns.distplot(target1['y'],kde=True,bins=50, color="blue", ax=ax[0,1])
sns.distplot(target1['width'],kde=True,bins=50, color="green", ax=ax[1,0])
sns.distplot(target1['height'],kde=True,bins=50, color="magenta", ax=ax[1,1])
locs, labels = plt.xticks()
plt.tick_params(axis='both', which='major', labelsize=12)
plt.show()
We can plot also the center of the rectangles points in the plane x0y. The centers of the rectangles are the points: $x_{c} = x + \frac{width}{2}$ and $y_{c} = y + \frac{height}{2}$.
We will show a sample of center points superposed with the corresponding sample of the rectangles.
from matplotlib.patches import Rectangle
fig, ax = plt.subplots(1,1,figsize=(7,7))
target_sample = target1.sample(2000)
target_sample['xc'] = target_sample['x'] + target_sample['width'] / 2
target_sample['yc'] = target_sample['y'] + target_sample['height'] / 2
plt.title("Centers of Lung Opacity rectangles (brown) over rectangles (yellow)\nSample size: 2000")
target_sample.plot.scatter(x='xc', y='yc', xlim=(0,1024), ylim=(0,1024), ax=ax, alpha=0.8, marker=".", color="brown")
for i, crt_sample in target_sample.iterrows():
ax.add_patch(Rectangle(xy=(crt_sample['x'], crt_sample['y']),
width=crt_sample['width'],height=crt_sample['height'],alpha=3.5e-3, color="yellow"))
plt.show()
tr = train_class_df[train_class_df['Target']==1]
centers = (tr.dropna(subset=['x'])
.assign(center_x=tr.x + tr.width / 2, center_y=tr.y + tr.height / 2))
ax = sns.jointplot("center_x", "center_y", data=centers, height=9, alpha=0.1)
_ = ax.fig.suptitle("Where is Pneumonia located?", y=1.01)
from sklearn.mixture import GaussianMixture
clf = GaussianMixture(n_components=2)
clf.fit(centers[['center_x', 'center_y']])
center_probs = clf.predict_proba(centers[['center_x', 'center_y']])
Z = -clf.score_samples(centers[['center_x', 'center_y']])
outliers = centers.iloc[Z > 17]
fig, ax = plt.subplots(1,1,figsize=(10,8))
centers.plot.scatter('center_x', 'center_y', c=Z, alpha=0.5, cmap='viridis', ax=ax)
outliers.plot.scatter('center_x', 'center_y', c='red', marker='x', s=100, ax=ax)
_ = ax.set_title('Where are the outliers?', fontsize=18)
image_sample_path = os.listdir('stage_2_train_images')[:5]
print(image_sample_path)
image_train_path = os.listdir('stage_2_train_images')
image_test_path = os.listdir('stage_2_test_images')
print("Number of images in train set:", len(image_train_path),"\nNumber of images in test set:", len(image_test_path))
print("Unique patientId in train_class_df: ", train_class_df['patientId'].nunique())
We confirmed that the number of unique patientsId are equal with the number of DICOM images in the train set.
Let's see what entries are duplicated. We want to check how are these distributed accross classes and Target value.
tmp = train_class_df.groupby(['patientId','Target', 'class'])['patientId'].count()
df = pd.DataFrame(data={'Exams': tmp.values}, index=tmp.index).reset_index()
tmp = df.groupby(['Exams','Target','class']).count()
df2 = pd.DataFrame(data=tmp.values, index=tmp.index).reset_index()
df2.columns = ['Exams', 'Target','Class', 'Entries']
df2
fig, ax = plt.subplots(nrows=1,figsize=(12,6))
sns.barplot(ax=ax,x = 'Target', y='Entries', hue='Exams',data=df2, palette='Set2')
plt.title("Chest exams class and Target")
plt.show()
Let's now extract one image and process the DICOM information.
import pydicom as dcm
samplePatientID = list(train_class_df[:3].T.to_dict().values())[0]['patientId']
samplePatientID = samplePatientID+'.dcm'
dicom_file_path = os.path.join("stage_2_train_images/",samplePatientID)
dicom_file_dataset = dcm.read_file(dicom_file_path)
dicom_file_dataset
def show_dicom_images(data):
img_data = list(data.T.to_dict().values())
f, ax = plt.subplots(3,3, figsize=(16,18))
for i,data_row in enumerate(img_data):
patientImage = data_row['patientId']+'.dcm'
imagePath = os.path.join("stage_2_train_images/",patientImage)
data_row_img_data = dcm.read_file(imagePath)
modality = data_row_img_data.Modality
age = data_row_img_data.PatientAge
sex = data_row_img_data.PatientSex
data_row_img = dcm.dcmread(imagePath)
ax[i//3, i%3].imshow(data_row_img.pixel_array, cmap=plt.cm.bone)
ax[i//3, i%3].axis('off')
ax[i//3, i%3].set_title('ID: {}\nModality: {} Age: {} Sex: {} Target: {}\nClass: {}\nWindow: {}:{}:{}:{}'.format(
data_row['patientId'],
modality, age, sex, data_row['Target'], data_row['class'],
data_row['x'],data_row['y'],data_row['width'],data_row['height']))
plt.show()
show_dicom_images(train_class_df[train_class_df['Target']==1].sample(9))
def show_dicom_images_with_boxes(data):
img_data = list(data.T.to_dict().values())
f, ax = plt.subplots(3,3, figsize=(16,18))
for i,data_row in enumerate(img_data):
patientImage = data_row['patientId']+'.dcm'
imagePath = os.path.join("stage_2_train_images/",patientImage)
data_row_img_data = dcm.read_file(imagePath)
modality = data_row_img_data.Modality
age = data_row_img_data.PatientAge
sex = data_row_img_data.PatientSex
data_row_img = dcm.dcmread(imagePath)
ax[i//3, i%3].imshow(data_row_img.pixel_array, cmap=plt.cm.bone)
ax[i//3, i%3].axis('off')
ax[i//3, i%3].set_title('ID: {}\nModality: {} Age: {} Sex: {} Target: {}\nClass: {}'.format(
data_row['patientId'],modality, age, sex, data_row['Target'], data_row['class']))
rows = train_class_df[train_class_df['patientId']==data_row['patientId']]
box_data = list(rows.T.to_dict().values())
for j, row in enumerate(box_data):
ax[i//3, i%3].add_patch(Rectangle(xy=(row['x'], row['y']),
width=row['width'],height=row['height'],
color="yellow",alpha = 0.1))
plt.show()
show_dicom_images_with_boxes(train_class_df[train_class_df['Target']==1].sample(9))
show_dicom_images(train_class_df[train_class_df['Target']==0].sample(9))
# PATH = ''
# vars = ['Modality', 'PatientAge', 'PatientSex', 'BodyPartExamined', 'ViewPosition', 'ConversionType', 'Rows', 'Columns', 'PixelSpacing']
# def process_dicom_data(data_df, data_path):
# for var in vars:
# data_df[var] = None
# image_names = os.listdir(PATH+data_path)
# for i, img_name in tqdm_notebook(enumerate(image_names)):
# imagePath = os.path.join(PATH,data_path,img_name)
# data_row_img_data = dcm.read_file(imagePath)
# idx = (data_df['patientId']==data_row_img_data.PatientID)
# data_df.loc[idx,'Modality'] = data_row_img_data.Modality
# data_df.loc[idx,'PatientAge'] = pd.to_numeric(data_row_img_data.PatientAge)
# data_df.loc[idx,'PatientSex'] = data_row_img_data.PatientSex
# data_df.loc[idx,'BodyPartExamined'] = data_row_img_data.BodyPartExamined
# data_df.loc[idx,'ViewPosition'] = data_row_img_data.ViewPosition
# data_df.loc[idx,'ConversionType'] = data_row_img_data.ConversionType
# data_df.loc[idx,'Rows'] = data_row_img_data.Rows
# data_df.loc[idx,'Columns'] = data_row_img_data.Columns
# data_df.loc[idx,'PixelSpacing'] = str.format("{:4.3f}",data_row_img_data.PixelSpacing[0])
# process_dicom_data(train_class_df,'stage_2_train_images/')
# train_class_df.sample(5)
# train_class_df.name = 'train_class_aug_df'
# train_class_df.csv_path = 'train_class_aug_df.csv'
# train_class_df.to_csv(train_class_df.csv_path)
# test_class_df = pd.read_csv('stage_2_sample_submission.csv')
# test_class_df.sample(5)
# test_class_df = test_class_df.drop('PredictionString',1)
# process_dicom_data(test_class_df,'stage_2_test_images/')
# test_class_df.name = 'test_class_aug_df'
# test_class_df.csv_path = 'test_class_aug_df.csv'
# test_class_df.to_csv(test_class_df.csv_path)
train_class_df = pd.read_csv('train_class_aug_df.csv')
test_class_df = pd.read_csv('test_class_aug_df.csv')
print(train_class_df.shape)
print(test_class_df.shape)
Let's check how many modalities are used. Both train and test set are checked.
print("Modalities: train:",train_class_df['Modality'].unique(), "test:", test_class_df['Modality'].unique())
The meaning of this modality is CR - Computer Radiography.
Let's check if other body parts than 'CHEST' appears in the data.
print("Body Part Examined: train:",train_class_df['BodyPartExamined'].unique(), "test:", test_class_df['BodyPartExamined'].unique())
View Position is a radiographic view associated with the Patient Position. Let's check the View Positions distribution for the both datasets.
print("View Position: train:",train_class_df['ViewPosition'].unique(), "test:", test_class_df['ViewPosition'].unique())
Let's get into more details for the train dataset. First, let's check the distribution of PA and AP.
get_feature_distribution(train_class_df,'ViewPosition')
Both AP and PA body positions are present in the data. The meaning of these view positions are:
Let's check, for the training data presenting Lung Opacity, the distribution of the window for both View Positions. We create a function to represent the distribution of the window centers and windows.
def plot_window(data,color_point, color_window,text):
fig, ax = plt.subplots(1,1,figsize=(7,7))
plt.title("Centers of Lung Opacity rectangles over rectangles\n{}".format(text))
data.plot.scatter(x='xc', y='yc', xlim=(0,1024), ylim=(0,1024), ax=ax, alpha=0.8, marker=".", color=color_point)
for i, crt_sample in data.iterrows():
ax.add_patch(Rectangle(xy=(crt_sample['x'], crt_sample['y']),
width=crt_sample['width'],height=crt_sample['height'],alpha=3.5e-3, color=color_window))
plt.show()
We sample a subset of the train data with Target = 1. We calculate as well the center of the windows with Lung Opacity. We then select from this sample the data with the two view position, to plot the window distribution separately.
target1 = train_class_df[train_class_df['Target']==1]
target_sample = target1.sample(2000)
target_sample['xc'] = target_sample['x'] + target_sample['width'] / 2
target_sample['yc'] = target_sample['y'] + target_sample['height'] / 2
target_ap = target_sample[target_sample['ViewPosition']=='AP']
target_pa = target_sample[target_sample['ViewPosition']=='PA']
plot_window(target_ap,'green', 'yellow', 'Patient View Position: AP')
plot_window(target_pa,'blue', 'red', 'Patient View Position: PA')
Let's check the Conversion Type data.
print("Conversion Type: train:",train_class_df['ConversionType'].unique(), "test:", test_class_df['ConversionType'].unique())
Both train and test have only WSD Conversion Type Data. The meaning of this Conversion Type is WSD: Workstation.
print("Rows: train:",train_class_df['Rows'].unique(), "test:", test_class_df['Rows'].unique())
print("Columns: train:",train_class_df['Columns'].unique(), "test:", test_class_df['Columns'].unique())
Only {Rows:Columns} {1024:1024} are present in both train and test
Let's examine now the data for the Patient Age for the train set.
g = sns.FacetGrid(col='Target', hue='PatientSex',
data=train_class_df.drop_duplicates(subset=['patientId']),
height=6, palette=dict(F="red", M="blue"))
_ = g.map(sns.distplot, 'PatientAge', hist_kws={'alpha': 0.3}).add_legend()
_ = g.fig.suptitle("What is the age distribution by gender and target?", y=1.02, fontsize=20)
Our CNN architecture stacks 4 convolutional layers (each one followed by a ReLU layer), then a pooling layer, then another few convolutional layers (+ReLU), then another pooling layer, and so on. The image gets smaller and smaller as it progresses through the network, but it also typically gets deeper and deeper (i.e., with more feature maps) thanks to the convolutional layers. At the top of the stack, a regular feedforward neural network is added, composed of a few fully connected layers (+ReLUs), and the final layer outputs the prediction (e.g., a softmax layer that outputs estimated class probabilities).
import csv
import random
import pydicom
from skimage import io
from skimage import measure
from skimage.transform import resize
import tensorflow as tf
from tensorflow import keras
from matplotlib import pyplot as plt
import matplotlib.patches as patches
# empty dictionary
PATH = ''
pneumonia_locations = {}
# load table
with open(os.path.join(PATH+'stage_2_train_labels.csv'), mode='r') as infile:
# open reader
reader = csv.reader(infile)
# skip header
next(reader, None)
# loop through rows
for rows in reader:
# retrieve information
filename = rows[0]
location = rows[1:5]
pneumonia = rows[5]
# if row contains pneumonia add label to dictionary
# which contains a list of pneumonia locations per filename
if pneumonia == '1':
# convert string to float to int
location = [int(float(i)) for i in location]
# save pneumonia location in dictionary
if filename in pneumonia_locations:
pneumonia_locations[filename].append(location)
else:
pneumonia_locations[filename] = [location]
# load and shuffle filenames
folder = PATH+'stage_2_train_images'
filenames = os.listdir(folder)
random.shuffle(filenames)
# split into train and validation filenames
n_valid_samples = 2560
train_filenames = filenames[n_valid_samples:]
valid_filenames = filenames[:n_valid_samples]
print('n train samples', len(train_filenames))
print('n valid samples', len(valid_filenames))
n_train_samples = len(filenames) - n_valid_samples
print('Total train images:',len(filenames))
print('Images with pneumonia:', len(pneumonia_locations))
class generator(keras.utils.Sequence):
def __init__(self, folder, filenames, pneumonia_locations=None, batch_size=32, image_size=256, shuffle=True, augment=False, predict=False):
self.folder = folder
self.filenames = filenames
self.pneumonia_locations = pneumonia_locations
self.batch_size = batch_size
self.image_size = image_size
self.shuffle = shuffle
self.augment = augment
self.predict = predict
self.on_epoch_end()
def __load__(self, filename):
# load dicom file as numpy array
img = pydicom.dcmread(os.path.join(self.folder, filename)).pixel_array
# create empty mask
msk = np.zeros(img.shape)
# get filename without extension
filename = filename.split('.')[0]
# if image contains pneumonia
if filename in pneumonia_locations:
# loop through pneumonia
for location in pneumonia_locations[filename]:
# add 1's at the location of the pneumonia
x, y, w, h = location
msk[y:y+h, x:x+w] = 1
# if augment then horizontal flip half the time
if self.augment and random.random() > 0.5:
img = np.fliplr(img)
msk = np.fliplr(msk)
# resize both image and mask
img = resize(img, (self.image_size, self.image_size), mode='symmetric')
msk = resize(msk, (self.image_size, self.image_size), mode='symmetric') > 0.5
img = np.array(img).astype(np.float32)
msk = np.array(msk).astype(np.float32)
# add trailing channel dimension
img = np.expand_dims(img, -1)
msk = np.expand_dims(msk, -1)
return img, msk
def __loadpredict__(self, filename):
# load dicom file as numpy array
img = pydicom.dcmread(os.path.join(self.folder, filename)).pixel_array
# resize image
img = resize(img, (self.image_size, self.image_size), mode='symmetric')
# add trailing channel dimension
img = np.expand_dims(img, -1)
return img
def __getitem__(self, index):
# select batch
filenames = self.filenames[index*self.batch_size:(index+1)*self.batch_size]
# predict mode: return images and filenames
if self.predict:
# load files
imgs = [self.__loadpredict__(filename) for filename in filenames]
# create numpy batch
imgs = np.array(imgs)
return imgs, filenames
# train mode: return images and masks
else:
# load files
items = [self.__load__(filename) for filename in filenames]
# unzip images and masks
imgs, msks = zip(*items)
# create numpy batch
imgs = np.array(imgs)
msks = np.array(msks)
return imgs, msks
def on_epoch_end(self):
if self.shuffle:
random.shuffle(self.filenames)
def __len__(self):
if self.predict:
# return everything
return int(np.ceil(len(self.filenames) / self.batch_size))
else:
# return full batches only
return int(len(self.filenames) / self.batch_size)
The ResNet is surprisingly simple. It starts and ends exactly like GoogLeNet (except without a dropout layer), and in between is just a very deep stack of simple residual units. Each residual unit is composed of two convolutional layers (and no pooling layer), with Batch Normalization (BN) and ReLU activation, using 3 × 3 kernels and preserving spatial dimensions (stride 1, SAME padding).
def create_downsample(channels, inputs):
x = keras.layers.BatchNormalization(momentum=0.9)(inputs)
x = keras.layers.LeakyReLU(0)(x)
x = keras.layers.Conv2D(channels, 1, padding='same', use_bias=False)(x)
x = keras.layers.MaxPool2D(2)(x)
# Added start
#x = keras.layers.Conv2D(channels, 1, padding='same', use_bias=False)(x)
#x = keras.layers.MaxPool2D(2)(x)
# Added End
return x
def create_resblock(channels, inputs):
x = keras.layers.BatchNormalization(momentum=0.9)(inputs)
x = keras.layers.LeakyReLU(0)(x)
x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)
x = keras.layers.BatchNormalization(momentum=0.9)(x)
x = keras.layers.LeakyReLU(0)(x)
x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)
#Added Start
x = keras.layers.BatchNormalization(momentum=0.9)(x)
x = keras.layers.LeakyReLU(0)(x)
x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)
#Added End
addInput = x;
print("Add input shape:", addInput.shape)
print("Resnet block input shape:", inputs.shape)
resBlockOut = keras.layers.add([addInput, inputs])
print("Resnet block out shape:", resBlockOut.shape)
out = keras.layers.concatenate([resBlockOut, addInput], axis=3)
print("concat block out shape:", out.shape)
out = keras.layers.Conv2D(channels, 1, padding='same', use_bias=False)(out)
print("mixed block out shape:", out.shape)
return out
def create_network(input_size, channels, n_blocks=2, depth=4):
# input
inputs = keras.Input(shape=(input_size, input_size, 1))
x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(inputs)
# residual blocks
for d in range(depth):
channels = channels * 2
x = create_downsample(channels, x)
for b in range(n_blocks):
x = create_resblock(channels, x)
# output
x = keras.layers.BatchNormalization(momentum=0.9)(x)
x = keras.layers.LeakyReLU(0)(x)
x = keras.layers.Conv2D(1, 1, activation='sigmoid')(x)
outputs = keras.layers.UpSampling2D(2**depth)(x)
model = keras.Model(inputs=inputs, outputs=outputs)
return model
# define iou or jaccard loss function
def iou_loss(y_true, y_pred):
y_true = tf.reshape(y_true, [-1])
y_pred = tf.reshape(y_pred, [-1])
intersection = tf.reduce_sum(y_true * y_pred)
score = (intersection + 1.) / (tf.reduce_sum(y_true) + tf.reduce_sum(y_pred) - intersection + 1.)
return 1 - score
# combine bce loss and iou loss
def iou_bce_loss(y_true, y_pred):
return 0.5 * keras.losses.binary_crossentropy(y_true, y_pred) + 0.5 * iou_loss(y_true, y_pred)
# mean iou as a metric
def mean_iou(y_true, y_pred):
y_pred = tf.round(y_pred)
intersect = tf.reduce_sum(y_true * y_pred, axis=[1, 2, 3])
union = tf.reduce_sum(y_true, axis=[1, 2, 3]) + tf.reduce_sum(y_pred, axis=[1, 2, 3])
smooth = tf.ones(tf.shape(intersect))
return tf.reduce_mean((intersect + smooth) / (union - intersect + smooth))
folder = PATH+'stage_2_train_images'
train_gen = generator(folder, train_filenames[0:1000], pneumonia_locations, batch_size=16, image_size=128, shuffle=False, augment=True, predict=False)
valid_gen = generator(folder, valid_filenames[1000:1200], pneumonia_locations, batch_size=16, image_size=128, shuffle=False, predict=False)
Gradient descent optimization algorithms like stochastic gradient descent (SGD) are most widely used to train neural networks. Recently, adaptive variants of gradient methods like Adagrad, Adam etc. have caught popularity and are being used by almost all the new architectures for training.Adam optimizer has become the de-facto standard for these algorithms.
Adam is an an algorithm for first-order gradient-based optimization of stochastic objective functions, based on adaptive estimates of lower-order moments. The method is straightforward to implement, is computationally efficient, has little memory requirements, is invariant to diagonal rescaling of the gradients, and is well suited for problems that are large in terms of data and/or parameters.
model = create_network(input_size=128, channels=16, n_blocks=2, depth=3)
model.compile(optimizer='adam',
loss=iou_bce_loss,
metrics=['accuracy', mean_iou])
print("model summary:", model.summary())
# cosine learning rate annealing
def cosine_annealing(x):
lr = 0.001
epochs = 25
return lr*(np.cos(np.pi*x/epochs)+1.)/2
learning_rate = tf.keras.callbacks.LearningRateScheduler(cosine_annealing)
# train_gen = generator(folder, train_filenames, pneumonia_locations, batch_size=16, image_size=128, shuffle=False, augment=True, predict=False)
# valid_gen = generator(folder, valid_filenames, pneumonia_locations, batch_size=16, image_size=128, shuffle=False, predict=False)
history = model.fit_generator(train_gen, validation_data=valid_gen, callbacks=[learning_rate], epochs=10, shuffle=True, verbose=2)
#history = model.fit_generator(train_gen, callbacks=[learning_rate], epochs=25, workers=4, use_multiprocessing=True)
plt.figure(figsize=(14,5))
plt.subplot(131)
plt.plot(history.epoch, history.history["loss"], label="Train loss")
plt.plot(history.epoch, history.history["val_loss"], label="Valid loss")
plt.legend()
plt.subplot(132)
plt.plot(history.epoch, history.history["accuracy"], label="Train accuracy")
plt.plot(history.epoch, history.history["val_accuracy"], label="Valid accuracy")
plt.legend()
plt.subplot(133)
plt.plot(history.epoch, history.history["mean_iou"], label="Train iou")
plt.plot(history.epoch, history.history["val_mean_iou"], label="Valid iou")
plt.legend()
plt.show()
model.save(r'Mixed_link_model.h5')
adm_accuracy = round(history.history.get('accuracy')[-1] * 100, 2)
adm_val_accuracy = round(history.history.get('val_accuracy')[-1] * 100, 2)
adm_loss = round(history.history.get('loss')[-1] * 100, 2)
adm_val_loss = round(history.history.get('val_loss')[-1] * 100, 2)
adm_mean_iou = round(history.history.get('mean_iou')[-1] * 100, 2)
adm_val_mean_iou = round(history.history.get('val_mean_iou')[-1] * 100, 2)
print(adm_accuracy, adm_val_accuracy, adm_loss, adm_val_loss, adm_mean_iou, adm_val_mean_iou)
for imgs, msks in valid_gen:
# predict batch of images
preds = model.predict(imgs)
# create figure
f, axarr = plt.subplots(2, 8, figsize=(20,8))
axarr = axarr.ravel()
axidx = 0
# loop through batch
for img, msk, pred in zip(imgs, msks, preds):
# plot image
axarr[axidx].imshow(img[:, :, 0])
# threshold true mask
comp = msk[:, :, 0] > 0.5
# apply connected components
comp = measure.label(comp)
# apply bounding boxes
predictionString = ''
for region in measure.regionprops(comp):
# retrieve x, y, height and width
y, x, y2, x2 = region.bbox
height = y2 - y
width = x2 - x
axarr[axidx].add_patch(patches.Rectangle((x,y),width,height,linewidth=2,edgecolor='b',facecolor='none'))
# threshold predicted mask
comp = pred[:, :, 0] > 0.5
# apply connected components
comp = measure.label(comp)
# apply bounding boxes
predictionString = ''
for region in measure.regionprops(comp):
# retrieve x, y, height and width
y, x, y2, x2 = region.bbox
height = y2 - y
width = x2 - x
axarr[axidx].add_patch(patches.Rectangle((x,y),width,height,linewidth=2,edgecolor='r',facecolor='none'))
axidx += 1
plt.show()
# only plot one batch
break
Stochastic Gradient Descent (SGD) is a variant of gradient descent in which the parameter update happens for every training example. Vanilla gradient descent computes the gradient of cost function for entire training dataset before performing one update and SGD gets rid of this redundancy by performing one update at a time.
Recent work shows that adaptive optimization techniques generalize much worse than stochastic gradient descent (SGD) even when the solutions have better training performance. Here we train the model with SGD optimizer and compare the model performance with Adam optimizer in the validation set.
model = create_network(input_size=128, channels=16, n_blocks=2, depth=3)
model.compile(optimizer='SGD',
loss=iou_bce_loss,
metrics=['accuracy', mean_iou])
print("model summary:", model.summary())
history = model.fit_generator(train_gen, validation_data=valid_gen, callbacks=[learning_rate], epochs=10, shuffle=True, verbose=2)
plt.figure(figsize=(14,5))
plt.subplot(131)
plt.plot(history.epoch, history.history["loss"], label="Train loss")
plt.plot(history.epoch, history.history["val_loss"], label="Valid loss")
plt.legend()
plt.subplot(132)
plt.plot(history.epoch, history.history["accuracy"], label="Train accuracy")
plt.plot(history.epoch, history.history["val_accuracy"], label="Valid accuracy")
plt.legend()
plt.subplot(133)
plt.plot(history.epoch, history.history["mean_iou"], label="Train iou")
plt.plot(history.epoch, history.history["val_mean_iou"], label="Valid iou")
plt.legend()
plt.show()
sgd_accuracy = round(history.history.get('accuracy')[-1] * 100, 2)
sgd_val_accuracy = round(history.history.get('val_accuracy')[-1] * 100, 2)
sgd_loss = round(history.history.get('loss')[-1] * 100, 2)
sgd_val_loss = round(history.history.get('val_loss')[-1] * 100, 2)
sgd_mean_iou = round(history.history.get('mean_iou')[-1] * 100, 2)
sgd_val_mean_iou = round(history.history.get('val_mean_iou')[-1] * 100, 2)
print(sgd_accuracy, sgd_val_accuracy, sgd_loss, sgd_val_loss, sgd_mean_iou, sgd_val_mean_iou)
result = pd.DataFrame(
{
'Optimizer': ['ADAM', 'SGD '],
'accuracy': [adm_accuracy, sgd_accuracy],
'val_accuracy': [adm_val_accuracy, sgd_val_accuracy],
'loss': [adm_loss, sgd_loss],
'val_loss': [adm_val_loss, sgd_val_loss],
'mean_iou': [adm_mean_iou, sgd_mean_iou],
'val_mean_iou': [adm_val_mean_iou, sgd_val_mean_iou]
})
result.sort_values(by=['accuracy', 'val_accuracy'], axis=0, ascending=True,
inplace=False, kind='quicksort', na_position='last')
result
How does your solution affect the problem in the domain or business?
Ans: Our solution will automatically locate lung opacities on chest radiographs.
What recommendations would you make, and with what level of confidence? Ans: Our solution with recomend the bounding boxses around the lung opacities on chest radiographs with 96.64% accuracy.